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A  LOCAL  INTERPOLATOR  DERIVED  FROM  THE 
DISCRETE  FOURIER  TRANSFORM 


1.  INTRODUCTION 

Interpolating  between  measured  data  points  is  a  common  signal  processing  problem.  Standard 
interpolation  methods  [1],  which  use  polynomial  fits  to  the  sampled  data,  perform  well,  but  even  better 
performance  is  desirable  for  critical  applications.  In  a  frame-differencing  signal  processor,  for  example, 
the  difference  between  two  successive  pictures  is  taken  to  "zero  out"  stationary  backgrounds  and  render 
faint  moving  targets  detectable.  Since  the  samples  from  the  two  pictures  (usually)  do  not  coincide,  values 
must  be  interpolated  between  samples  of  the  second  picture  to  compare  to  the  first.  The  more  accurately 
this  can  be  done,  the  less  clutter  will  occur  in  the  differenced  image. 

In  the  ideal  case,  when  there  is  a  large  number  of  samples  and  the  Nyquist  sampling  criterion  is 
satisfied,  the  Nyquist  reconstruction  formula  (NRF)  supplies  the  answer  [2,  e  pluribus  unum].  But  the 
theoretical  NRF  uses  an  infinite  sum,  which  must  be  truncated  in  some  manner  for  practical  application 
and  generally  performs  poorly  near  the  edge  of  a  picture  (for  example),  where  the  truncation  becomes 
severe. 


Polynomial  spline  techniques  lack  computational  simplicity  because  they  require  a  global  solution: 
interpolated  values  at  any  point  depend  on  data  points  that  are  far  away  [3].  Furthermore,  if  a 
higher-order  spline  is  desired,  a  completely  different  solution  must  be  applied. 

The  usual  approach  to  purely  local  inteipolation  is  to  apply  a  low-order  polynomial  fit  to  the  data 
in  a  short,  sliding  window  consisting  of  an  even  number  of  points  (e.g.,  N  =  2,  4,  or  6)  to  obtain 
interpolated  values  between  the  two  central  points.  Figure  1  illustrates  the  problem  for  N  =  4.  Since  only 
a  few  points  are  used,  interpolations  can  be  done  very  close  to  the  edge  of  the  picture.  The  discrete 
Fourier  transform  (DFT)  interpolator  described  here  is  also  local  in  nature,  hence  is  computationally  easy 
to  apply.  Furthermore,  it  uses  the  same  formula  for  any  number  of  data  points  and,  depending  on  the  data 
sampling  rate,  can  give  results  superior  to  the  local  polynomial  interpolators. 


+  SAMPLED  DATA.  )xn\ 
■  RESAMPLING  POINT 


I 
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| 

*4*- 
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I 

i  o  i  j  .i  ; 

SAMPLE  NUMBER 

i  _  _ 

Fig.  1  —  The  solid  curve  represents  the  analog  signal.  An 
interpolated  value  at  the  shifted  position  s  must  be  found  from  the 
sampled  data  and  compared  to  the  analog  signal 
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2.  THE  DFT  INTERPOLATOR 

Appendix  A  shows  how  the  DFT  and  the  Fourier  shift  theorem  (FST)  can  be  used  to  motivate 
a  definition  for  an  analog  of  the  NRF  for  a  limited  number  of  data  points.  In  one  dimension,  for  an  even 
number  of  points,  the  result  is  Eq.  (A7),  which  is  always  applied  with  n  =  M 2-1;  as  with  the  local 
polynomials,  interpolations  are  done  only  between  the  two  central  points  in  a  sliding  window,  as  indicated 
in  Fig.  1.  Putting  this  restriction  into  Eq.  (A7)  and  appropriately  relabeling  the  summation  index  to 
conform  to  Fig.  1,  we  have  an  interpolation  formula  to  apply  to  a  set  of  points  {;cn}, -M2 -1)  <n<  M2, 
to  obtain  an  interpolated  value  for  the  point  xs,  0  <  s  <  1 : 

jc  =  1  Y  x  sin  [*(*  ~  ”)]  (1) 

*  A  („=_W2-1)  "  tan  [tt(s  -  n)IN]’ 

when  N  is  even.  The  forms  for  odd  N  and  for  two  dimensions  are  given  by  Eqs.  (A8)  and  (A9).  Note 
that  since  N  'tdnlir(s -n)/N]  -*  n(s-n )  as  N  -*  oo,  Eq.  (1)  approaches  the  NRF  in  this  limit. 

Equation  (1)  has  the  form  of  a  convolutional  interpolation  formula 

Nf2 

xs  =  xn  /4„(s), 

n=-(Nt. 2-1) 


where  the  interpolation  coefficients  A„(s)  are  obtained  from  an  interpolation  function  F(x)  by 

An(s)  =  F(n-s).  (3) 

The  interpolation  function  for  DFT-A  ( N  even)  is  obviously 


F(x)  = 


sin(7r.r) 

N  tan(TTJc/iV) 


w. 


(4) 


where  RN  is  the  rectangle  function  of  width  N  centered  at  the  origin;  that  is, 
Rfjx)  =  1  for  | jc  |  <  M2,  R^/jc)  =  0  otherwise.  Equation  (4),  as  applied  in  Eq.  (1),  now  stands  as  a 
candidate  interpolation  function  to  be  compared  to  conventional  polynomial  interpolation  functions. 

Previous  applications  of  the  DFT  to  the  interpolation  problem  have  centered  on  a  "zero  padding" 
method  in  frequency  space  due  to  Schafer  and  Rabiner  [4]  and  recently  investigated  by  Fraser  [5].  This 
method  is  appropriate  if  a  number  of  additional  samples  evenly  spaced  between  existing  samples  is 
desired.  An  example  would  be  a  requirement  for  ten-times-denser  resampling  with  nine  equally  spaced 
new  samples  between  each  pair  of  existing  samples.  But  the  technique  is  not  well  suited  to  generating  a 
single  new  sample  at  an  arbitrary  location.  The  approach  presented  here  is  well  suited  to  this  problem 
and  can  be  shown  to  give  the  same  interpolated  value  as  the  zero  padding  technique  by  manipulating 
Fraser’s  Eq.  (5).  This  can  more  readily  be  seen  from  his  Eq.  (6),  which  is  derived  from  his  Eq.  (5),  but 
his  Eq.  (6)  unfortunately  contains  an  error  that  prevents  the  calculated  interpolated  value  from  being  a 
pure  real  number  for  nonintegral  interpolation  points  (whereas  his  Eq.  (5)  does  guarantee  a  real  value). 
Once  this  error  is  corrected  (see  Appendix  B  of  this  report),  Fraser’s  Eq.  (6)  is  identical  (in  content)  to 
this  report’s  Eq.  (A5). 
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Finally,  an  alternate  form  of  the  interpolation  function  that  does  not  run  the  computational  risk 
of  a  zero  denominator  is 


m 


2 

N 


N/Z-l 

+  cos(7rx)  +  2  cos(2 nirx/N) 

fl-\ 


for  N  even,  and 


(5) 


F( x)  =  I 
N 


+ 


(N- 1)/2 


2E 

n  -  1 


cos  (2  mrx/N) 


R^x) 


(6) 


for  N  odd. 

Combining  Eqs.  (5)  or  (6)  with  Eqs.  (2),  (3),  and  (4)  shows  that  the  DFT  interpolators  are  a 
special  case  of  trigonometric  interpolation  [6] .  They  are  equivalent  to  using  the  lowest-order  trigonometric 
polynomial  that  can  fit  N  sample  points  over  the  interval  (-M2,M2)  when  the  sample  points  are  at 
-N/2  4-  1,  -M2  +  2,  ...,M 2,  for  N  even  and  at  -(N  -  l)/2,  -(N  -  l)/2  +  1,  ...,(A  -  l)/2  for  N  odd. 

3.  TESTING  THE  INTERPOLATOR 

The  coefficients  of  an  Mpoint  DFT  are,  of  course,  the  coefficients  of  the  first  M2  4-  1  terms 
(including  DC,  see  Section  4)  in  a  Fourier  series  describing  the  data.  If  the  data  actually  consisted  of  just 
those  discrete  frequencies  used  by  the  DFT,  then  the  interpolation  would  be  exact.  To  see  that  the 
interpolation  formula  works  well  on  more  realistic  signals,  even  for  small  N  (N  =  4,6,8),  it  will  be 
compared  to  polynomial  interpolators  applied  to  a  representative  point  spread  function  (PSF). 

A  PSF  is  used  as  a  test  waveform  because  i‘  is  the  most  rapidly  varying  signal  that  will  occur  in 
most  scenes  and,  therefore,  the  hardest  to  interpolate.  (An  example  of  a  scene  that  can  produce  a  more 
rapidly  varying  signal  is  a  positive-contrast  point  source  located  very  close  to  a  negative-contrast  point 
source,  but  such  occurrences  are  rare.)  The  representative  PSF  chosen  is  a  unit-height  Gaussian,  with 
standard  deviation  a  given  in  units  of  the  sample  spacing.  Thus,  by  changing  a,  we  can  see  how  well  an 
interpolator  works  as  a  function  of  sampling  rate. 

The  DFT  interpolators  are  obtained  from  Eq.  (1)  with  the  appropriate  value  of  N  and  will  be 
compared  to  polynomials  that  use  the  same  number  of  data  points  to  perform  an  interpolation.  Thus 
DFT-4,  DFT-6,  and  DFT-8  will  be  compared  to  third-,  fifth-,  and  seventh-order  polynomials, 
respectively.  The  third  order  (cubic)  interpolator  used  is  that  described  by  Wolberg  [7]  and  extensively 
investigated  by  Park  and  Schowengerdt  [8],  with  the  parameter  <x  chosen  to  be  -0.5.  This 
"image-independent  optimal  choice"  [8]  emphasizes  high  fidelity  at  low  frequencies.  The  fifth-order 
(quintic)  interpolator  was  derived  by  Schaum  [9],  using  a  more  general  method.  The  seventh-order 
interpolator  will  be  standard  Lagrange  interpolation  [  1],  denoted  LF-8  (Lagrange  function  for  N  =  8). 
The  formulas  for  these  interpolators  are  given  in  Appendix  C.  Figure  2  plots  the  interpolation  functions 
for  all  these  interpolators. 
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Fig.  2  —  The  DFT  interpolation  functions  from  Eq.  (4)  for  N  = 
4,  6,  8.  The  polynomial  interpolation  functions  obtained  from 
Appendix  3  and  Eq.  (3)  are  included  for  comparison. 


Figure  3  shows  the  results  for  a  =  1  and  worst-phased  sampling  (samples  miss  the  peak  of  the 
curve  by  the  maximum  amount).  Figure  4  shows  best-phased  sampling  (a  sample  occurs  at  the  peak  of 
the  curve).  The  DFT  interpolators  fit  die  peak  of  a  Gaussian  better  than  the  polynomials  but  tend  to  ring 
more  at  the  foot  of  the  curve.  Errors  calculated  over  the  intervals  shown  in  Figures  3  and  4  are  given  for 
various  a’s  in  Table  1. 

The  following  remarks  about  Table  1  emphasize  worst-phased  sampling  because  it  is  bound  to 
occur  at  many  places  in  real  data.  For  a  =  1.5  (1.5  samples  per  standard  deviation),  the  polynomials  are 
superior  to  their  DFT  counterparts,  but  the  differences  are  small;  with  well-sampled  data,  all  interpolators 
perform  well.  For  a  =  0.5,  none  give  acceptable  results  (though  the  polynomials  may  be  marginally 
adequate  for  best-phased  sampling).  The  advantage  of  the  DFT  interpolators  becomes  apparent  for  the 
intermediate  case  of  a  =  1.  DFT-4,  6,  and  8  give  somewhat  better  results  than  the  corresponding 
polynomials  for  root-mean-square  (rms)  errors  and  substantially  better  results  for  maximum  errors. 
(DFT- 10,  not  shown,  yields  a  small  further  improvement,  but  the  point  of  diminishing  returns  has  been 
reached.) 

The  greater  superiority  of  the  DFT  interpolators  with  respect  to  maximum  errors  than  with 
respect  to  rms  errors  (as  shown  in  Table  1)  can  be  important.  In  the  frame-differencing  signal  processor 
alluded  to  earlier,  for  example,  reducing  maximum  errors  will  more  quickly  reduce  the  number  of  false 
alarms  for  a  given  threshold  setting  than  reducing  the  rms  error  will. 

The  fact  that  DFT-7  gives  worse  results  than  DFT-6  is  dealt  with  in  Section  4.  It  is  worth  noting 
that  if  the  NRF  is  truncated  to  ten  points  or  fewer  and  normalized  to  have  unit  response  at  zero 
frequency,  its  performance  is  substantially  worse  than  any  of  the  interpolators  shown  here.  It  is  also 
worth  noting  that  cr  =  1  is  a  somewhat  critical  value.  If  the  PSF  is  sampled  substantially  more  often  than 
once  per  standard  deviation,  the  choice  of  interpolator  becomes  unimportant.  If  it  is  sampled  substantially 
less  often,  none  perform  well. 
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c  i 
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SAMPLE  NUMBER 


Fig.  3  —  The  solid  curves  arc  unit-height  Gaussians  with  standard 
deviation  equal  to  the  sample  spacing.  The  dashed  lines  show  how 
well  different  interpolators  perform  in  reproducing  the  Gaussians 
from  the  marked  samples,  with  worst-phased  sampling. 
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Fig.  4  —  Same  as  Fig.  1,  but  with  best-phased  sampling.  In  this 
case,  all  the  interpolators  except  cubic  fit  the  Gaussian  peak  almost 
perfectly. 
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Table  1  —  Interpolation  Errors* 


<T  = 

Interpolation 

Method 

Worst  Phase 

Best  Phase 

max 

rms 

max 

rms 

0.5 

Cubic 

0.33 

0.084 

0.022 

0.006 

Quintic 

0.30 

0.079 

0.043 

0.014 

LF-8 

0.28 

0.076 

0.071 

0.022 

DFT-4 

0.28 

0.077 

0.066 

0.021 

DFT-6 

0.26 

0.070 

0.091 

0.034 

DFT-7 

0.23 

0.064 

0.13 

0.063 

DFT-8 

0.25 

0.067 

0.11 

0.041 

1 

Cubic 

0.050 

0.013 

0.032 

0.012 

Quintic 

0.030 

0.008 

0.020 

0.007 

LF-8 

0.022 

0.006 

0.013 

0.005 

■HI9 

Ml— 

0.C11 

0.029 

0.010 

DFT-6 

0.016 

0.006 

0.014 

0.005 

DFT-7 

0.084 

0.026 

0.073 

0.026 

DFT-8 

0.010 

0.004 

0.008 

0.003 

1.5 

Cubic 

0.013 

0.004 

0.012 

0.C04 

Quintic 

0.005 

0.002 

0.006 

0.002 

LF-8 

0.003 

0.001 

0.003 

0.001 

DFT-4 

0.021 

0.009 

0.021 

0.009 

DFT-6 

0.010 

0.004 

0.009 

0.004 

DFT-7 

0.076 

0.030 

0.072 

0.030 

DFT-8 

0.006 

0.002 

0.006 

0.002 

*Entries  are  maximum  and  root-mean-square  errors  over  the  intervals  shown  in  Figs.  1 
and  2  for  Gaussians  with  the  indicated  standard  deviations 


6 


NRL  REPORT  9392 


Testing  an  interpolator  on  real  data  is  very  desirable  but  seldom  possible;  an  interpolated  value 
can  be  calculated,  but  the  true  value  to  which  it  should  be  compared  is  usually  not  available.  To  show 
that  the  superior  performance  of  the  DFT  interpolators  is  not  limited  to  Gaussian  functions,  they  were 
compared  to  the  quintic  polynomial  interpolator  on  accurately  simulated  data.  This  was  done  by  using 
a  scene  simulator  that  produces  images  of  the  Earth  as  seen  from  a  satellite,  with  50-  to  100-m  resolution. 
The  PSF  for  the  simulated  scene  is  close  to,  though  not  exactly,  a  Gaussian  (but,  of  course,  very  little 
of  the  scene  produces  as  hard-to-fit  a  signal  as  the  PSF,  so  we  should  expect  good  interpolation 
performance  for  a  <  1).  On  scenes  with  little  structure,  all  the  interpolators  performed  well.  On  a  highly 
structured  scene  (including  an  aerial  view  of  a  city)  with  a  *  0.5,  all  interpolators  gave  about  equally 
poor  results.  With  a  =  0.75,  all  gave  reasonably  good  results,  but  the  DFT  interpolators  performed 
better  than  the  polynomials  by  about  the  margins  shown  in  Table  1  for  a  =  1  with  worst-phased 
sampling.  With  a  >  1,  all  the  interpolators  performed  well;  quintic  was  slightly  better  than  DFT-6,  but 
not  quite  as  good  as  DFT-8. 

The  lesson  of  the  simulated  imagery  and  of  Table  1  is  that  the  DFT  interpolators  give  better 
performance  at  lower  sampling  rates  (i.e.,  smaller  a)  than  the  polynomials. 

4.  DISCUSSION 

To  see  why  the  DFT-based  interpolation  formulas  work  well,  we  first  consider  in  more  detail  the 
fact  that  interpolations  are  done  only  between  the  two  center  points  of  a  sliding  window. 

Interpolation  with  the  DFT  and  FST  is  sometimes  done  globally  by  transforming  the  entire  data 
string,  multiplying  the  Fourier  coefficients  by  the  appropriate  phase  factors  (as  is  done  with  the  gk  in 
Appendix  1),  and  transforming  back  to  real  space  to  obtain  a  set  of  interpolated  samples  [10,11],  This 
has  the  advantage  of  doing  the  interpolations  between  all  data  points  simultaneously  (as  opposed  to 
between  the  two  central  points  only)  but  often  gives  disappointing  results  because  of  edge  discontinuities. 

Discontinuities  at  the  edges  of  a  data  set  are  a  common  problem  with  Fourier  techniques;  if  the 
data  (or  its  derivatives')  do  not  wrap  around  smoothly,  a  Fourier  analysis,  which  assumes  cyclically 
repeating  data,  "sees"  a  discontinuity  at  the  edge.  Near  any  discontinuity,  a  Fourier  series  converges 
slowly  and  imperfectly,  an  effect  known  as  Gibbs’  phenomenon  and  covered  in  standard  texts  [12],  This 
problem  is  generally  dealt  with  by  various  "windowing"  methods  (e.g.,  von  Hann  [  i  3]  or  Hamming  [14]). 
But  these  methods  are  imperfect  (they  impose  their  own  structure  on  the  image  and  are  not  band-limited) 
ai  )  impractical  in  a  small  data  set.  Figure  5,  in  which  six-  and  seven-point  DFT’s  are  applied  to  a  ramp, 
shows  the  effect  of  an  edge  discontinuity  when  interpolations  are  done  between  all  the  data  points 
simultaneously.  The  six-point  DFT  can  do  accurate  interpolations  between  the  two  central  points;  in  this 
interval,  the  effect  of  the  edge  discontinuity  is  quite  small.  But  interpolations  for  the  seven-point  DFT 
are  strongly  affected  in  all  intervals.  From  this  example,  we  expect  a  DFT  interpolator  to  give  much 
better  results  for  an  even  number  of  points  than  for  an  odd  number,  and  the  comparison  of  DFT-7  with 
DFT-6  in  Table  1  shows  this  to  be  the  case. 

Once  a  DFT-based  interpolation  formula  is  considered  to  be  a  standard  interpolation  process 
(applying  only  to  the  center  interval  of  a  sliding  window),  it  becomes  apparent  that  the  main  reason  for 
its  success  is  its  good  frequency  response. 

An  N-point  DFT  describes  its  input  data  in  terms  of  a  few  "allowed  frequencies,"  i.e.,  frequencies 
for  which  an  integral  number  of  cycles  fit  exactly  into  the  interval  covered  by  N  sample  points.  The 
maximum  frequency  contained  in  the  DFT  is  fm  =  1/2,  (i.e.,  one  cycle  per  two  samples),  regardless  of 
the  number  of  samples.  Thus  even  a  few-point  DFT  directly  measures  the  highest  frequency  available  in 
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Fig.  5  —  The  dashed  curves  show  the  values  interpolated  between 
samples  of  a  ramp  by  six-point  (circles)  and  seven-point  (squares) 
DFT  interpolators.  The  former  reproduces  the  central  interval  in  the 
ramp  very  well,  while  the  latter  docs  not  work  well  in  any  interval. 


the  sampled  scene  (the  Nyquist  frequency)  as  well  as  the  lowest  (the  DC  component).  Of  course, 
the  few-point  DFT  has  very  poor  frequency  resolution;  the  difference  between  successive  frequencies 
is  A f  =  1  IN,  so  there  are  only  M2  nonzero  frequencies  in  the  range  from  zero  frequency  to  the  Nyquist 
frequency.  The  DFT  approximates  intermediate  frequencies  (such  as  one  half-way  between  two  of  the 
allowed  frequencies)  by  a  sum  of  the  allowed  frequencies. 

The  above  remarks  are  for  even  N ,  which  is  the  case  of  interest  here.  For  the  sake  of 
completeness,  the  corresponding  quantities  for  odd  N  are  Af  -  1/M  as  before,  but  fm  =  (1-1  /AO/2; 
therefore,  the  number  of  nonzero  frequencies  is  ( N  -  l)/2.  Thus  a  six-point  and  a  seven-point  DFT  both 
contain  three  nonzero  frequencies;  the  difference  is  that  the  six-point  (or  any  even-point)  DFT  contains 
only  the  cosine  phase  of  its  highest  frequency,  while  the  seven-point  (or  any  odd-point)  DFT  also  contains 
the  sine  phase,  as  shown  in  Fig.  6. 


'.5  *  tir-  tj  *r  * }  *r*t»  *rt_r  *■  > 


Fig.  6  —  If  three  complete  cycles  of  a  sinusoid  are  covered  by  six 
samples  (circles),  only  the  cosine  phase  (solid  line)  yields  nonzero 
values,  the  sine  phase  (dashed  line)  goes  undetected.  For  seven 
samples  (squares),  both  signals  are  detected.  The  "missing"point  on 
the  right  belongs  to  the  next  six-  or  seven-point  interval. 


But  even  though  a  DFT  interpolator  has  poor  frequency  resolution,  it  interpolates  its  allowed 
frequencies  perfectly,  and  these  frequencies  are  evenly  distributed  across  the  spectrum.  This  results  in 
a  more  uniform  frequency  response  than  is  possible  with  the  polynomials  that  can  be  optimized  to  give 
better  Fits  at  some  frequencies  than  others  [8]  but  cannot  reproduce  any  nonzero  frequency  exactly  (since 
their  outputs  can  only  be  a  piece-wise  continuous  polynomial,  not  a  true  sinusoid).  We  will  follow 


8 


NRL  REPORT  9392 


standard  practice  in  taking  the  Fourier  transform  of  the  interpolation  function  as  a  useful  measure  of  the 
interpolator’s  frequency  response  (however,  see  below  for  further  remarks  on  this  point).  The  Fourier 
transform  of  F(x)  (which  is  easily  found  from  Eq.  (5))  is 

<  N/2-1 

G{k)  =  l{sinc  [7V(A:-l/2)]  +  sine  {A/(*+l/2)]}  ♦  £  sine  [N(k-nlN)\,  (7) 

7  n— (N/2-1) 

for  even  N,  where  sinc(x)  *  sin (irx)l(icx). 

Figure  7  shows  the  Fourier  transforms  of  the  cubic,  quintic,  LF-8,  and  DFT-4,  6,  and  8 
interpolation  functions  (all  the  transforms  are  real  and  symmetric,  so  only  the  positive  halves  are  shown). 
Note  that  the  DPT  interpolators  have  unit  response  at  their  allowed  frequencies  (i.e.,  unity  at/  =  1/4  for 
DFT-4,  at  1/6  and  1/3  for  DFT-6,  1/8,  1/4,  and  3/8  for  DFT-8),  except  at/  =  1/2,  where  the  response 
of  0.5  reflects  the  fact  that  only  the  cosine  phase  of  this  frequency  is  sampled.  The  maximum  value  at 
lower  frequencies  is  1.03.  Also  note  that  the  DPT'  interpolators  have  broader-band  response  than  the 
comparable  polynomials.  The  ideal  frequency  response  is  the  rectangle  function  /?,  (unity  between  +  1/2, 
zero  otherwise;  this  is  the  Fourier  transform  of  the  NRF’s  sine  interpolation  function).  The  DFT 
interpolation  functions  approach  this  limit  as  their  order  increaces.  The  approach  is  especially  rapid  at 
low  frequencies,  which  is  important  because  the  power  spectral  density  of  virtually  all  imagery  has  a 
much  higher  value  at  low  frequencies  than  at  high  frequencies. 

Interpreting  the  Fourier  transform  of  the  interpolation  function  as  the  frequency  response  of  the 
interpolator  is  not  as  clear-cut  as  one  would  like.  First,  the  term  "frequency  response”  must  be  employed 
carefully  for  an  interpolator  because  it  is  not  well  defined;  an  input  sinusoid  produces  an  output  that  is 
not  a  sinusoid.  This  is  true  for  the  polynomials  at  any  nonzero  frequency  and  for  the  DFTs  at  any 
frequency  other  than  the  allowed  frequencies.  Second,  as  shown  in  Appendix  D,  if  an  interpolator 
reproduces  an  input  frequency  exactly,  then  the  Fourier  transform  will  have  the  value  unity  at  that 
frequency,  but  the  converse  is  not  true.  For  example,  the  parameters  of  the  cubic  interpolation  function 
can  be  chosen  so  that  its  Fourier  transform  has  unit  value  at  some  point  in  addition  to  zero  [8],  but  since 
the  result  of  an  interpolation  is  a  piece-wise  cubic  function,  an  input  sinusoid  at  that  frequency  is  not 
reproduced  exactly. 


L4T-)  urn  dftb 
I  „  1.  _  t  l-O' 

*  f'v.- 

C  Q  X 

\ 


0  11  1  .i  i  I 

I  I  'v  "i 

Fig.  7  —  The  Fourier  transforms  of  the  cubic 
(C),  quintic  (Q),  LF-8,  and  DFT  interpolation 
functions.  The  Nyquist  frequency  is/  =  1/2. 
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5.  CONCLUSION 

DFT-based  interpolators,  when  applied  locally  to  the  center  of  a  sliding  window,  have  been 
shown  to  be  a  valuable  addition  to  the  signal  processor’s  repertoire.  On  highly  structured  data,  they  can 
extend  good  interpolation  performance  to  lower  data  sampling  rates  than  is  possible  with  local  polynomial 
interpolators. 
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APPENDIX  A 


Derivation  of  Interpolation  Formula 

The  form  of  the  Nyquist  reconstruction  formula  for  a  limited  number  of  data  points  is  exhibited 
below  (Eqs.  (A7),  (A8),  or  (A9)).  A  sketch  of  the  derivation  is  given  for  a  one-dimensional  discrete 
Fourier  transform  (DFT)  for  even  N.  The  changes  in  the  derivation  for  odd  N  are  then  indicated,  and  the 
generalization  to  two  dimensions  is  presented. 

For  a  set  of  N  equally  spaced  samples  {*„},  0  <,n<,  AM  from  a  continuous  function  fiu)  (i.e., 
xn  ~  An))>  the  coefficients  of  the  corresponding  DFT  are  (following  the  conventions  of  Oppenheim  and 
Schafer  [Al]) 


N- 1 


Thus, 


gk  =  Y  x„  e-fitkn/N  Q^k^N-l 

n= 0 


v-iEi. 

"  *= 0 


(Al) 


(A2) 


Recognizing  that  gNk  =  g\for  k  ^  1  (and  therefore,  that  for  even  N,  gN/2  =  g*N/2  =  real  number),  we 
can  write  this  as 


x„  =  — 
"  N 


JV/2-1 


go  +  W-0"  +  2 Re  £  gk 


(A3) 


The  Fourier  shift  theorem  (FST)  [A2]  states  that  if  F(k)  is  the  (continuous)  Fourier  transform  of 
J{x),  then  the  Fourier  transform  of  fix  +  s)  is  F(k)&xp(iks).  The  interpolation  problem  consists  of  finding 
a  value  for  fin  +  s )  given  {xn} .  When  the  DFT  is  expressed  in  the  form  Eq.  (A3),  comparison  of  the 
FST  and  the  interpolation  problem  suggests  the  following  definition  of  a  new  set  of  DFT  coefficients: 

8k  =  8k  ej2*ks,N  for  k  <.  NI2-1, 
and 


8'm  =  8m  Re(.ejws)  =  gm  cos(xs), 


(A4) 


where  s,  0  £  s  ^  1  denotes  the  shift.  The  definition  of  g  'n/2  is  chosen  to  assure  that  the  coefficient  of 
the  highest  frequency  in  the  DFT  is  real,  as  it  must  be  for  N  even,  if  it  is  to  be  the  DFT  of  a  real-valued 
(as  opposed  to  complex-valued)  sequence  of  data.  By  defining  {xn4J}  to  be  the  sequence  of  real  numbers 
that  will  result  from  the  inverse  DFT  of  the  coefficients  defined  in  Eq.  (A4)  and  substituting  their  values 
in  terms  of  the  gk,  we  have  a  candidate  expression  for  an  interpolated  sample: 
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=  — 


1_ 

N 


Nfl-l 


So 


gm(- 1)"  cos(7T5)  +  2 Re  £  gk  eW"'VN 


(A5) 


Whether  xn+s  closely  approximates  ./(n  +  s)  remains  to  be  seen  (for  n  =  M2  - 1,  it  does,  see  text).  We 
can  now  substitute  for  the  g' s  from  (Eq.  (Al)),  use  the  facts  that  (-iy  m  cos(ir^)  =  cos[ir(n  +  s  -  m)\  and 


m- i 


y  uk  =  u  x-.-u _ =  u 

h  i  - « 


m-i 


-  ,.NI 4 


u(N-2)/4  _  U-(N- 2)/4 

u1/2-«-1/2  ’ 


(A6) 


with  u  =  exp[/2ir(«  +  5  -  m)/A]  to  obtain,  for  even  A, 

.  /v-i 


sin  [7r(n+5-/n)] 


±  y  x  _ _ 

N  ~z o  m  tan[ir(n+s-/n)/A) 


=  _ (7r~y)  52  cot  [ir(n*s-m)IN].  (A7) 

A  m=0 

For  odd  A,  M 2-1  is  replaced  by  (A  -  l)/2  in  Eq.  (A3)  and  g0  is  the  only  term  that  appears  outside  the 
summation.  The  rest  of  the  derivation  is  unchanged.  The  result  for  odd  A  is 

jcntJ  =  y  xm(-l rm  CSC  Mn+s-m)/N\.  (A8) 

A  m=0 

The  generalization  to  two  dimensions  is  obviously 


xn  +  j,m + r 


sin(7r.v)  sin(irr) 

W 


AM  M- 1 

EE  xk/-w 

k=Q  /= 0 


,n*m-k-l 


cot  [it(n+s-k)IN]  cot  [7r(m  +r-[)IM] .(A9) 


Here  the  one-dimensional  interpolation  is  applied  sequentially  in  each  direction.  If  either  (both)  of  A  and 
M  are  odd,  then  the  corresponding  cotangent  function(s)  is  (are)  replaced  by  cosecant(s). 
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APPENDIX  B 


Correction  of  Earlier  Work 

This  appendix  corrects  an  error  in  Fraser’s  1989  paper  on  fast  Fourier  transform  interpolation 
[Bl];  it  must  be  read  in  conjunction  with  that  paper  to  be  comprehensible. 

N  will  be  taken  to  be  even  (for  JVodd,  Fraser’s  paper  contains  no  error),  so  Fraser’s  N'l 2  can 
simply  be  replaced  by  N/2.  To  see  that  Fraser’s  Eq.  (6)  contains  an  error,  the  special  case  U(N/2)  =  real 
number,  all  other  U(k)  =  0,  will  be  considered.  Then  Fraser’s  Eq.  (5),  which  correctly  gives  the  values 
interpolated  from  the  original  sequence,  which  is  denoted  by  u(i),  i  =  0,1,2 -  1,  becomes 

v(i')  =  V(M2)exp \j2iri ' (NI2)/L]  *  V(L-M2)expl/27ri'(L-M2)/L] 

=  0.5  *  U(NI2)  [exp(jri'IM)  +  exp(-yViVM)] 


=  U(NI2)  cos(ir/7M),  i'  =  0,1,2,... £-1, 


(Bl) 


where  Fraser’s  Eqs.  (3)  and  (4)  and  L  =  M-N have  been  used.  In  Eq.  (Bl),  the  interpolated  values  are 
always  real.  However,  for  the  same  special  case,  Fraser’s  Eq.  (6),  which  is  ostensibly  obtained  from  his 
Eq.  (5)  by  letting  i'/Af  -*x  (that  is,  letting  i'/M  become  a  continuous  variable),  gives 

v(Mx)  =  U(N/2)  exp  (jirx) 

=  U(NI2)exp(jiri' IM)  for  x  =  i'lM  (B2) 

which  is  not  real  for  most  i\  Thus,  there  must  be  an  error  in  Fraser’s  Eq.  (6).  After  putting  x  =  i'lM 
in  Fraser’s  Eq.  (5),  a  few  straight-forward  steps  lead  to 

m-l  AV2-I 

v(Mx)  =  ^  U(k)exp{j2irkx/N)  +  J2  U(k)*  Qxp(-j2irkx/N) 

Jt=0  k=\ 

+  0.5f/(M2)[exp(/7rx)  +  exp  (-jirx)), 

where  use  has  been  made  of  the  fact  that  U(NI 2)  must  be  real  because  it  is  the  N/2  coefficient  of  an 
/V-term  DFT  (with  N  even)  of  real  data.  This  is  the  correct  version  of  Fraser’s  Eq.  (6).  It  can  be  seen  to 
be  identical  in  content  (though  not  in  the  1  IN  normalization  convention)  to  this  report’s  Eq.  (A5)  by 
setting  x  =  n  +  s. 
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APPENDIX  C 


Polynomial  Interpolation  Formulas 


The  formulas  used  for  the  polynomial  interpolators  are  given  below. 


Cubic 


Ouintic 


=  i  E  **  An> 
L  n«-l 


t  =  -L  y'  xn  b 

24  i  " 

^  n*-2 


where 


A_ j  =  s(l  -  s )2 
A0  =  2-s2(5-3s) 

Ax  =  s(l  +  4s  -  3s2) 
A2  =  -  s2(l  -s) 


fl_2  =  s(2-s)-s3(9  -  13^  +  5s2) 

=  -  16s(l  -s)  +  s3(39  -  645  +  25 52) 
fi0  =  24  -  30s2  -  s3(70  -  126s  +  50s2) 
5,  =  16s(l  +  s)  +  s3  (66  -  124s  +  50s2) 
B2  =  -  s(2+s)  -  s3(33  -  61s  +  25s2) 

B3  =  s3(7  -  12s  +  51) 


Standard  Lagrange  interpolation  [Cl]  consists  of  fitting  an  (N  -  1)  order  polynomial  through  N  data 
points.  Specialized  to  the  problem  of  interest  here,  this  yields  LF  -  8: 


4  4 

"  E  *n  Cn>  With  C„  =  J] 

it* -3  m*-3 

m#n 


(s-m) 

(»-»)’ 
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APPENDIX  D 


Interpolation  Errors  for  Pure  Sinusoids 


First  we  define  two  error  terms  Ee  and  E,  by  the  equations 

Nfl 

52  cos(2xJbn)  F(s-m)  =  cos(2xfo)  +  Ec(k,s) 

m— (M/2-1) 

and 


(Dl) 


N/ 2 


£ 

m— (M/2-1) 


sin(2 rkm)  F(s-m)  =  sin(2xto) 


Es(k,s), 


(D2) 


where  0  <  s  <  1.  The  left-hand  sides  of  Eqs.  (Dl)  and  (D2),  where  F(x)  is  the  interpolation  function 
from  Eq.  (3),  will  be  recognized  from  Eq.  (1)  as  the  interpolated  values  of  cosine  and  sine  inputs, 
respectively.  Thus,  Ec  and  Es  are  the  errors  resulting  from  the  interpolation.  In  what  follows,  the  first 
equality  is  the  Fourier  transform  of  the  interpolation  function  (which  is  an  even  function  [e.g.,  Eq.  (4)], 
so  its  Fourier  transform  has  no  imaginary  component).  The  second  equality  is  obtained  by  the  change  of 
variables  x  -*  m  -  s  for  each  value  of  m.  The  third  equality  is  then  found  by  expanding  cos[2x£(.v  ~m)] , 
using  F(m  -  s)  -  F(s  -  m),  and  substituting  from  Eqs.  (Dl)  and  (D2). 

G{k)  =  J  "  F( x)  cos  (2rkx)dx  (D3> 


i  m 

[  J2  cos[2x^(s-m)]  F(s-m)ds 

■»0  m— (M/2-1) 


(D4) 


=  1  f  1  (Ec(k,s)  cos(2xfa)  +  Es(k,s)  sin(2 -xks))ds.  (E>5) 


This  is  the  desired  result:  if  Ec  -  Es  =  0  for  some  k  (e.g.,  for  a  DFT  interpolator  at  an  allowed 
frequency),  then  G(k)  =  1.  But  G(k)  can  be  unity  with  Ec,  Es  *  0,  if  the  integral  in  Eq.  (D5)  is  zero. 
There  is,  unfortunately,  no  simple  relation  between  the  error  of  the  interpolator  and  the  nonunity  of  the 
Fourier  transform  of  the  interpolation  function. 
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